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OPTIMAL  ESTIMATION  OF  MEASUREMENT  BIAS 


INTRODUCTION .  A  Kalman  filtering  program  has  been  developed  by  the 
authors  to  provide  a  Best  Estimate  of  Trajectory  (BET)  for  flight  tests 
conducted  at  White  Sands  Missile  Range  (WSMR).  This  optimal  filtering 
program  combines  measurements  from  radar,  fixed  camera,  cinetheodolite , 
dovap,  velocimcter,  and  accelerometer  which  are  optimally  weighted  using 
on-line  estimates  cf  the  measurement  variances.  One  of  the  most 
important  considerations  in  developing  a  BET  technique  is  to  account  for 
inconsistencies  produced  by  bias  errors  in  the  measurements  As  a  matter 
of  fact,  a  BET  is  not  very  useful  unless  these  measurement  errors  have 
been  accounted  for.  It  was  for  this  reason  that  this  research  project 
began.  As  a  result  of  this  project,  we  now  have  a  very  effective  tech¬ 
nique  included  in  the  BET  program  to  estimate  the  measurement  bias  errors. 

Since  our  BET  program  uses  a  Kalman  filter  it  was  desirable  that  the 
bias  estimation  technique  be  developed  within  the  framework  of  the 
Kalman  filter  theory.  There  is  a  natural  way  of  including  bias  terms 
in  the  Kalman  filter.  One  merely  adds  an  additional  state  variable  for 
each  bias  term  to  be  considered  and  forms  the  optimal  estimate  of  the 
biases  in  the  same  way  as  for  the  other. states  of  the  system.  This 
technique  is  fine  for  cases  where  there  are  only  a  few  bias  terms  to  be 
considered.  However,  a  typical  application  of  our  PET  program  has  a 
large  number  of  measuring  instruments  involved.  For  example^  a  LANCE 
flight  test  might  have  two  radars,  28  dovap  receivers,  eight  fixed 
cameras,  and  eight  cinetheodolites.  Considering  only  one  bias  error 
per  measurement  this  results  in  66  additional  state  variables  to  be 
estimated.  If  the  dimension  of  dynamic  state  is  nine,  we  would  then 
have  to  compute  filtered  estimates  for  7b  state  variables.  An  ordinary 


Kalman  filter  program  using  a  75  dimensional  state  vector  is  computa¬ 
tionally  p.ohibitive  at  the  present  time.  Thus,  we  must  develop  some 
technique,  other  than  the  straightforward  method  of  augmenting  the 
filter  state  vector,  for  estimating  measurement  bias  errors  within  the 
Kalman  filter  framework.  The  research  problem  may  be  stated  concisely 
as  follows:  Develop  a  computationally  feasible  method  within  the  Kalman 
filter  framework  for  estimating  the  biases  of  all  measuring  instruments 
participating  in  a  WSMR  mission. 

A  preliminary  search  of  the  literature  reveals  two  possible  tech¬ 
niques  which  might  be  applicable  to  the  bias  estimation  problem.  A 
third  technique  developed  by  the  author  was  also  considered.  Finally, 
during  the  course  of  the  research  an  additional  technique  was  found. 

Of  the  four  methods  considered,  all  of  which  were  computationally 
feasible,  three  were  discarded  either  because  of  numerical  difficulties 
or  excessive  errors  in  the  bias  estimates.  The  remaining  method,  which 
was  published  by  B.  Friedland  (in  Reference  (1))  was  developed  and 
extensively  evaluated  for  WSMR  applications.  Several  extensions  of 
Friedland' s  method  were  made  so  that  it  would  be  applicable  to  the 
instrument  bias  problem  at  WSMR.  Also,  this  report  employs  a  derivation 
different  than  Friedlands, 

The  evaluation  of  this  instrument  bias  estima'.ion  technique  was 
performed  using  simulated  measurement  data  from  a  nearly  ballistic 
trajectory.  The  simulated  instrumentation  included  four  fixed  cameras, 
six  cinetheodolites,  one  accelerometer,  and  two  radars  with  doppler. 
Thus,  there  were  a  total  of  twenty-nine  bias  states  estimated  in 
addition  to  nine  dynamic  states  of  the  basic  filter.  The  technique  was 
also  tested  using  real  data  from  a  similar  trajectory  having  the  same 
instrumentation.  The  simulated  data  has  bias  and  noise  added  to  each  of 
the  exact  measurements.  Evaluation  of  time  varying  as  well  as  constant 
bias  were  considered  in  the  evaluation. 
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BIAS  ESTIMATION  WITH  THE  KALMAN  FILTER.  We  will  consider  only  the 
discrete  case  of  the  Kalman  filter.  The  following  equations  descri?  3 
the  discrete  Kalman  filter  as  used  in  this  report. 

DYNAMIC  STATE  EQUATION.  The  model  of  the  process  we  are  observing  is 
represented  by  the  linear  difference  equation 

x(k+l)  =  A(k)  x(k)  +  u^Ck)  (1) 

(n*l)  .  (n*n)  (n*l)  (n*l) 

where  is  the  state  vector  of  the  process  and  u^Ck)  is  a  random  vector 
representing  our  uncertainties  in  how  well  the  model  defined  by  the 
homogeneous  portion  of  the  difference  equation  acually  represents  the 
system.  We  assume  that 


and 

E[\(k)"*U)]  = 

OBSERVATION  EQUATION.  At  discrete  instants  of  time  t^  we  have  vector 
valued  observations  z(i)  available.  The  z(i)  are  assumed  to  be  linearly 
related  to  the  state  by 

2<i)  =  H(i)  x(i)  +  G(i)  b(i)  +  v(i)  (2) 

(mxi)  (m*n)  (n*l)  (m*p)  (p*l)  (m*l) 

where  v(i)  is  a  random  vector  of  measurement  errors  with  mean  zero  and 
covariance 


E  jv£i)vT(  j)J  =  R(i)6ij 


Sg»-Ng 


■1  ■ 
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b(k)  is  a  vector  of  measurement  biases  which  obey 


b(k+l)  =  B(k)b(k)  +  u^k) 


where. u.  (k)  is  a  zero  mean  random  vector  with  covariance 

D 


*  vk)5u 


(3) 


By  adjoining  the  vector  b  to  the  vector  x  we  form  the  augmented  state 
vector  y 


y  = 


(4) 


Then  the  dynamics  of  y  are 


y<k+X>  *  r(k)v(U  *  u(k) 


(5) 


where 


F(k)  = 


A(k) 

n*n 


B(k) 

pxpj 


and 


u(k )  = 


\<k> 


ub(k) 


^  r-.  ^  ,  |  aaassaBiaMfeaaiitaaa&bi  a  rrrhWrWirr  m  - at  aaaaflalhig^ttiiaBiSiau 


where 


The  posterior  covariance  P(k+1)  may  be  defined  as 

P(k+1)  =  jV^k+lJk)  +  LT(k+l)R"1(k+l)L(k+l)J"1  (12) 

Direct  implementation  on  a  digital  computer  of  the  augmented  Kalman 
filter  equations  defined  by  (7)  thru  (12)  is  computationally  prohibitive 
when  the  dimension  of  the  bias  vector  b  is  large.  However,  by  utilizing 
suitable  restrictions  on  the  bias  dynamics,  Priedland  (Reference  (1))  was 
able  to  decouple  the  augmented  Kalman  filter  so  that  a  dynamic  state 
x*(k)  and  a  bias  b(k)  are  separately  estimated  and  then  the  optimal 
state  estimate  x(k)  is  computed  by 

x(k)  =  x*(k)  +  T(k)b(k)  (13) 


where  T(k)  is  an  nxp  matrix.  The  state  estimate  x*  is  computed  by 
assuming  all  biases  to  be  zero.  F;  will  find  that  there  are  certain 
restrictions  which  must  be  placed  on  the  form  of  the  augmented  filter 
in  order  that  the  decoupling  be  possible.  One  such  restriction  will  be 
that  Q^,  the  covariance  of  the  stochastic  term  in  the  bias  equation,  be 
zero.  The  development  is  not  restricted  to  measurement  biases;  biases 
in  assumed  constants  in  the  state  dynamics  may  also  be  included. 


Priedland  approached  the  decoupling  of  the  Kalman  filter  by  trans¬ 
formation  of  the  discrete  Ricatti  equation.  The  derivation  given  in  this 
report  will  be  based  upon  an  examination  of  the  conditions  for  which 
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the  decoupling  of  the  state  estimates  specified  by  (13)  is  possible. 

FILTER  DECOMPOSITION.  A  general  form  for  a  discrete  time,  linear 
recursive  filter  for  the  augmented  state  vector  may  he  written  in  the 
form  of  (9) 

y(k)  =  y(k|k-l)  +  W(k)^z(k)-H(k)x(k|k-l)-G(k)b(k|k-l)j  (14) 

where  for  the  present  we  will  consider  W(k)  to  be  an  arbitrary  gi'.in 
matrix.  The  augmented  state  estimate  y(k)  may  be  decomposed  into 

A  A 

estimates  x(k)  and  b(k) 

x(k)  =  x(k|k-l)  +  Wx(k)^z(k)-H(k)x(k|k-l)-G(k)b(k|k-l)j  (15) 

•x(kjk-l)  =  A(k-l)x(k-l) 

b(k)  s  b(k|k-l)  +  Wb(k)^z(k)>H(k)x(k|k-l)-G(k)b(k|k-l)J  (16) 

A  A 

b(k|k-l)  =  B(k-l)b(k-l) 


where  W(k)  in  (14)  is  . 


W(k)  = 


Vk> 


Vk) 


How  let  x*(k)  be  the  state  estimate  that  would  be  obtained  if  all 
biases  are  assumed  to  be  zero.  We  will  call  this  the  zero-bias 
estimate.  The  recursive  estimation  equations  for  x*(k)  are 
x*(k)  -  x*(k| k-1)  +  W, (k) (z(k)  -  H(k)x*(k|  k-1)) 

(1/ 

x*(k|  k-1)  -  A(k-l)xMk-l) 
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How  we  ask  the  question:  Under  what  conditions  can  we  decompose 

A 

the  optimal  estimate  x(k)  as 

x(k)  =  x*(k)  +  T(k)b(k) 


Denote  the  residual  by  r(k|k-l). 

r(k|k-l)  =  z(k)  -  H(k)x(k|k-1)  -  G.(k)b(k|k-1)  (1€ 

If  the  decomposition  (13)  holds,  substitution  of  (13)  with  k  replaced 

A 

by  k-1  into  (18)  results  in  an  alternative  expression  for  r(k|k-l) 

r(kjk-l)  =  z(k)  -  H(k)x*(k|k-1)  -  S(k)b(k-  -  )  (IS 


where 


S(k)  s  G(k)B(k-l)  +  H(k)A(k-l)T(k“l) 


How  substituting  (16),  (17),  (19),  and  (20)  into  (13)  we  find 

x(k)  =  x(kjk-l)  +  |w*(k)  +  T(k)Wb(k)Jr(k|k-l) 


+  T(k)B(k)  +  Wx(k)S(k)  -  A(k-l)T(k-l)  b(k-l) 


Note,  if^we  assume  that  (15)  and  (16)  define  the  optimal  linear 
filter,  r(kjk-l)  and  b(k-l)  must  be  linearly  independent  random  vectors. 
In  this  case  the  following  conditions  must  hold  in  (21), 


yk>  =  wx(k)  +  T(k)wb(k) 


T(k)  =  A(k-l)T(k-l)B‘1(k-l)  -  H*(k)S(k)B“1(k-l) 


8 


By  examination  of  (21)  we  easily  see  that  an  arbitrary  linear  filter  is 
also  decomposable  into  the  form  specified  by  (13),  if  the  above  condi¬ 
tions  hold.  However,  these  conditions  are  not  necessary  for  decomposi¬ 
tion  of  a  general  linear  filter.  In  addition  to  the  conditions  specified 
by  (22),  which  we  will  call  the  gain  condition  and  (23),  the  recursion 
equation  for  T(k),  we  must  also  have 

x(0)  =  :t“(0)  +  T(0)b(0)  (24) 

since  we  assume  that  the  decomposition  holds  for  all  k.  A  natural 

choice  for  the  initial  conditions  to  satisfy  (24)  will  in  many  cases  be 
*  * 

x(0)  «  x  (0),  T(0)  *  0,  i.e.  we  assume  that  the  biases  have  no  initial 
effect  on  the  optimal  estimate. 

The  conditions  expressed  by  (22),  (23),  and  (24)  are  the  basic 
requirements  for  a  general  linear  filter  to  be  decomposed. 


RESTRICTIONS  IMPOSED  BY  THE  GAIN  CONDITION.  The  consequences  of  the 
gain  condition  stated  in  (22)  will  be  examined  when  both  the  estimators 

ft 

for  x.  and  x(k),  b(k)  are  Kalman  filters.  The  gain  for  the  Kalman 

K  it 

filter  giving  the  zero-bias  estimate  x  can  be  written  as 


W*(k)  =  P*(k)HT(k)R~1(k)  (25) 

ft 

The  matrix  P  (k)  in  this  case  (with  biases  actually  present)  may  be 
defined  as 


P*(k)  =  El/  x*(k)-x(k)+T(k)b(k)j 


(k)-x(k)+T(k)b(k) 


We  will  now  examine  the  relations  between  the  various  Kalman  covariance 
matrices  Px(k),  Px(k),  and  P^k).  Rewrite  the  decomposition  equation 

as 
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x  (k)  -  x(k)  +  T(k)b(k)  x(k)-x(k)  +  T(k)  (b(k)-b(k))j  (27) 

Then  multiplying  this  equation  by  its  transpose  and  taking  expected 
values  we  find 


Px<k>  =  Vk>  +  T(k)Pb(k)TT(k)  -  P^OO^Ck)  -  T(k)Pbx(k)  (28) 


where 


Pxb(k)  =  Ef<^k>-x(k)>^(k)-b(k))T),  Pbx(k)  =  P^(k) 

A  similar  result  for  P  (k|k-l) 

P*(k|  k-1)  *  Px(k| k-1)  +  A(k-l)T(k-l)Pb(k-l)T7(k"l)>^_1) 


-pxb(k!k-i)TT(k-i)AWA(k-i)T(k-1)pbX(k-1) 


The  gain  for  the  augmented  Kalman  filter  is 


Vk)  Px(k)  Pxb(k>l|H(k> 


Vk)J  [Pbx(k)  Pb<k>J  LG  (k>. 


R”1(k) 


Substituting  the  gain  condition  into  the  first  component  of  (31) 


fpx(k)HT(k)  +  rxb(k)GT(k)yi(k)  =  H*(k)  +  T(k)Wfa(k) 


Then  substituting  the  second  component  of  (31)  for  W.  (k)  and  (25)  for 
*  15 
Wx(k),  we  find  the  satisfaction  of  the  gain  condition  requires 


Pxb(k)(GCk)+H(k)T(k)]T  =  T(k)Vk)(G<k)+H<k)T(k))1 


Thus,  the  gain  condition  holds  if  we  impose  the  requirement  that 


Pxb(k)  =  T<k)Pbfk) 


We  have  translatedthe  gain  condition  into  a  covariance  restriction,  but 
we  must  continue  to  examine  the  covariance  relations  to  determine  the 
meting  of  (34).  The  relation  between  gain  and  covariance  in  the 
filter  can  be  written  as 


P(k)  »  | 

Jl-W(k)L(k)Jp(k|k-l) 

• 

px(k) 

Pxb(k)l  [l“Vk)H(k) 

-Wx(k)G(k)“ 

Vk) 

Pb(k)j  ”  -Wb(k)H(k) 

I-W.  (k)G(k) 
b  J 

Px(k!k-1) 

p^tkik-i) 

P^lclk-l) 

Pb(k|k-1) 

(35) 


(36) 


Then 

Pxb(k)  =  (I“Vk)H<k))Pxb(k^k"1)  “  Vk)G(k)Pb(Mk~1)  (37) 

Pb(k)  =  ^I-Wb(k)G(k))pb(k|k-l)  -  Wb(k)H(k)Pxb(k|k-l)  (38) 

How  suppose  that  the  covariance  condition  (34)  holds  at  k-1,  P^k-l)  = 
T(k-l)  Pb(k-1) 
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Then 


P^kjk-l)  =  A(k-l)T(k-l)Pb(k-l)DT(k-l)  (39) 

Pb(k|k-1)  =  B(k-l)Pb(k-l)BT(k-l)  +  Qb(k-1)  (40) 

% 

Substituting  (39)  and  (40)  in  (37)  and  (30)  and  using  (20)  and  (22)  gives 


Prt(k)  =  T(k) ^  S(k-l)-Wb(k)S(k)j  Pb(k-l)BT(k-l) 

-  [\f*(k)+T(k)Wb(k)j  G(k)Qb(k>l) 

Pb(k)  -  ^  B(k-l)-Wb(k)S(k) j  Pb(k-l)BT(k-l) 

+  (  i-wb(k)s(k)  )Qb(k-l) 


(41) 


(42) 


Clearly,  the  condition 

Pxb(k)  =  T(k)pb(k) 


holds  if 

Qb(k-1)  =  0 


Thus,  we  cannot  have  a  stochastic  term  in  the  dynamics  of  the  bias  states. 
There  does  not  seem  to  be  any  way  of  removing  this  restriction.  This 
restriction  may  lead  to  a  divergence  problem  in  the  bias  estimation,  but 
fortunately  we  can  treat  this  problem  in  another  way  which  will  be 
discussed  later. 
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SUMMARY.  The  following  formulas  summarize  the  decomposition  of  the 
optimal  linear  filter  for  the  process  described  by  (1)  thru  (6). 


ZERO-BIAS  FILTER 

i 

•  X*(kJ  k-1)  =  A(k-l)x*(k-l) 

P*(k|k-1)  =  A(k-l)P*(k-l)AT(k-l)  +  Q  (k-1) 

**  X 


P*(k) 

X 

S  (  P^klk-l)  +  H(k)R"1(k)HT(k)  j"1 

x*(k) 

=  x*(kjk-l)  +  W*(k) ^  Z(k)-H(k)x*(k|k-1) j 

W*(k) 

X 

=  P*(k)HT(k)R"l(k) 

BIAS  riLTER 

S(k) 

=  G(k)B(k-l)  +  H(k)T(klk-l) 

C(k) 

=  H(k)P*(kj  k-l)HT(k)  +R(k) 

(43) 

b(k  k-1) 

A 

=  B(k-l)b(k-l) 

Pb(k|k-1) 

=  B(k-l)Pb(k-l)BT(k-l) 

(44) 

Pb(k) 

=  ^^(kjk-l)  +ST(k)C"1(k)S(k)J"1 

(45) 

b(k) 

=  b(kjk-l)  +Wfa(k)  ^Z(k)-H(k)x*(k|k-l)-S(k)b(k|k-l)J 

wb(k) 

=  Pb(k)ST(k)c"1(k) 

(46) 
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OPTIMAL  STATE  ESTIMATE 


T(kjk-l)  =  A(k-l)T(k-l) 

T(k)  =  T(k|k-1)  -  W*(k)S(k) 

x(k)  =  x*(k)  +  T(k)b(k) 

px(k)  =  P*(k)  +  T(k)pb(k)TT(k)  (47) 

The  numbered  equations  in  the  summary  have  not  been  derived  but  are  easily 
obtainable.  Equation  (44)  is  merely  (40)  with  Qb(k)=0,  and  (47)  comes 
from  (28)  with  pxb(k)=T(k)pb(k) .  Equations  (43,  (45),  and  (46)  can  be 
derived  by  substituting  the  relations  between  covariances  in  (12). 

EXTENSION  TO  THE  NONLINEAR  CASE.  The  extension  of  the  bias  esti¬ 
mation  procedure  to  the  case  where  the  state  dynamics  are  nonlinear  and 
the  observational  equations  are  nonlinear  functions  of  the  stale  car.  be 
performed  in  any  one  of  the  ways  by  which  the  Kalman  filter  is  usually 
extended  to  the  nonlinear  case.  Let  the  discrete  time  nonlinear 
dynamics  be  represented  by 

x(k)  =  f(x(k-l),k-l)  +  ux(k)  (48) 

and  let  the  nonlinear  observation  equations  be  written  as 

z(k)  =  h(x(k),k)  +  G(x(k),k)b(k)  +  v(k)  (49) 

also  we  will  employ  the  additional  notation 

B<k)  *  [slkr(h<x(k,'k),G(x<k)’k))l^(k|k-i).b(k|ic.i) 
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G(k)  =  G(x(k|k-l),k) 


A(k-l)  =  f  -al^— -I’K"1!]  A 

L  9*k-l  -I  x(k-'l) 


Since  the  extended  Kalman  filter  is  not  optimal  we  view  the  decomposition 
as  applying  to  an  arbitrary  linear  filter.  Also,  the  decomposition  of 
the  extended  filter  is  not  an  exact  procedure  as  in  the  linear  case  but 
holds  only  approximately  since  an  additional  linearization  is  required 
to  derive  the  gain  condition  and  recursive  definition  of  T(k). 


.The  extended  Kalman  filter  is 


x(k)  =  x(k|k-lHWx(k)(*(k).h(x(k|k-l),l)-G(x(k|k-l),k) 
b(k|k-l) j 

b(k)  =  b(k|k-J)+Wb(k>!  z(k)-h^x(k|k-l),k  J  -C^x(k|k-l),kj 

b(k|k-l) j 


where 


x(kjk-l)  =  f ^  x(k-l),k-l j 
b(kjk-l)  =  B  (k-l)b(k-l) 


\<» 

Vk> 

P  (k)  P  ,  (k) 
x  xb 


Pbx(k)  Pb(k) 


HT(k) 


GT(k) 


R_1(k) 
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x*(k|k-l)  =  f(x*(k-l),k) 

W*(k)  =  P*(k)HT<f(k)R*’1(k) 

X  X 

« 

P*(k)  =  [^P*'1(k}k-l)+H*T(k)R_1(k)H*(k)]  “-1 

P*(k|k-1)  =  A*(k-l)P"’(k-l)A*r(.k-l)  +  Q  (k-1) 

X  XX 

where 


(51) 

(52) 

(53) 

(54) 


The  derivation  of  the  conditions  under  which  the  decomposition  approxi¬ 
mately  holds  proceeds  exactly  as  in  the  linear  case.  The  additional 
linearizations 

f(x(k-l))=f(x*(k-l),k)  +  A*(k-l)T(k-l)b(k-l)  (55) 

h(x(k|k-l),k)=h(x“(k|k-l),k)  +  HS‘(k)A*(k-l)T(k~l)b{k~l)  (56) 

G(x(k|k-l),k)-G(x*(k|k-l),k)  =  G*(k)  (57) 


are  required  in  obvious  places.  The  details  of  the  derivation  will  not 
be  given  but  the  results  of  the  decomposition  are  summarized  below. 
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ZERO-BIAS  TILTER 


x*(k|k-l)  =  f(x*(k|k-l),k) 

P*(k|k-1)  =  A(k-l)P*(k-l)AT(k-l)  +  Qv(k-1) 

X  X 


F>) 

P*rl(k|k-1)  +  H*(k)R"1(k)H*T(k) j  _1 

x*(k) 

=  x*(kjk-l)  +  ^^(k)  (z(k)-h(x  (k|k-l),k)) 

Wx(k) 

=  P*(k)H*T(k)R_1(k) 

H*(k) 

-  ph(x,k)']  * 

[_  3x  J  x  (k|k-l) 

BIAS  FILTER 

• 

S*(k) 

s  G*(k)B(k-l)  +  H*(k)T(k|k-l) 

(58) 

C*(k) 

=  H*(k)P*(kJk-l)H*T(k)  +  R(k) 

(59) 

b(klk-l)  =  B(k-l)b(k-l)  (60) 

Pb(k|k-1)  =  B(k-l)Pb(k-l)BT(k-l)  (61) 

Pb(k)  =  (P^CkJk-l)  +  S*T(k)C*'1(k)S*(k))";L  (62) 

b(k)  =  b(k|k-l)  v?fc(k)(z(k)-h(xK(k{k-l)-S1'(k)^(kJk>l))  (63) 

Wb(k)  =  Pb(k)S*T(k)c*"1(k)  (64) 

^  # 

•  * 

IS 


APPROXIMATELY 

OPTIMAL  STATE  ESTIMATE 

T(k| k-1) 

=  A(k-lWk-l) 

(65) 

T(k) 

=  T(k|k-1)  -  wY'(k) 

(66) 

i(k) 

=  (k)  +  T(k)b(k) 

(67) 

px<x) 

=  P*(k)  +  T(k)P.  (k)TT(k) 
x  b 

(68) 

BIAS  ESTIMATION  V7ITH  THE  FADING  MEMORY  KALMAN  FILTER.  The 
restriction  on  the  filter  decomposition  imposed  by  the  gain  condition 
that  there  be  no  state  noise  on  the  bias  state  variable,  i.e.  0^00=0 
may  lead  to  a  filter  divergence  problem  as  previously  noted.  Thus  in 
the  bias  estimation  technique  we  need  to  develop  some  other  method  for 
devaluing  the  effect  of  old  observations  on  the  bias  state  estimates. 

The  divergence  problem  arises  from  numerical  errors  in  computation  and 
from  errors  in  modeling  the  dynamics  and  covariance  matrices.  The 
divergence  caused  by  numerical  errors  is  most  easily  controlled  by 
employing  some  standard  numerical  analysis  techniques  and  reformulating 
the  Kalman  filter  in  terms  of  the  square  roots  of  the  covariance 
matrices,  see  Reference  (  2  ).  He  use  this  reformulation  and  also 
process  only  scalar  observations  in  our  BET  program  to  control  numerical 
errors.  The  treatment  of  mismodeling  errors  is  ■  jnsiderably  more 
difficult.  The  mismodeling  errors  may  either  be  intentional,  e.g., 
modeling  a  measurement  bias  as  a  single  state  variable  when  we  have  a 
much  more  realistic  bias  model  available  containing  several  state 
variables,  or  inadvertantly  when  we  do  not  have  a  realistic  model  avail¬ 
able.  Intentional  mismodeling  may  arise  because  we  consider  some  terms 
in  the  bias  model  to  be  small  and  we  do  not  wish  to  unnecessarily 
complicate  the  filter.  We  may  also  delete  some  more  important  terms 
from  the  bias  model  when  it  is  obvious  that  the  geometry  of  the 
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instrumentation  is  not  sufficient  to  estimate  these  terms. 


The  decomposition  procedure  will  not  permit  most  of  the  conventional 
methods  of  treating  errors  caused  by  mismodeling,  e.g,,  we  have  already 
seen  that  no  state  noise  on  the  bias  variables  is  permitted.  Also,  the 
method  of  directly  overweighting  the  observations  is  also  excluded  by  the 
decomposition.  Another  method  for  devaluing  the  effect  of  old  observations 
on  the  estimates  is  the  fading  memory  Kalman  filter  which  we  were  able  to 
show  is  allowed  by  the  decomposition.  The  fading  memory  filter  weights 
the  observations  exponentially  according  to  their  age.  It  was  first  used 
in  the  Kalman  framework  by  Fagin  {5]  and  more  recently  by  Tarn  and 
Zaborsky  [6],  and  Sorensen  and  Sachs  [7].  Sorensen  and  Sachs  showed  that 
the  fading  factors  need  not  be  constant  as  previously  used  and  also 
exhibited  some  other  previously  unreported  properties  of  the  fading  memory 
Kalman  filter. 

The  basic  idea  of  the  fading  memory  filter  is  that  observations  should 
be  assigned  increasingly  larger  variances  as  they  become  older.  Let  tn  be 
the  current  time  at  which  an  estimate  is  desired.  At  tn  we  model  the 
dynamics  and  the  observations  (using  the  augmented  state  vector  y)  by 


y(i,n)  »  A(i)y(i-1- ,n)  +  ux(i,n)  (69) 

z(i)  «*  L(i)y(i,n)  +  v(i,n)  (70) 


The  observation  noise  and  state  noise  covariances  are  defined  as 


n-1  \ 


1  / 

v(i,n)vT(i,n)  =  R(i,n)  »  R(i)exp  J  c 

L  J  V  j-1  V . 


(71) 
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E  P ux< i *P >ux< i *n ) 1  =  Qx^i,n*  =  Q(i)exp(  \  c 
*-  J  \  j=i+l 


(72) 


where  Q(i)  and  R(i)  are  the  usual  covariance  matrices.  The  scalars  c^>0 
are  called  the  fading  factors.  From  (71)  and  (72)  we  can  see  that  the 


covariances  assigned  to  the  state  noise  and  observation  noise  become 
larger  as  their  time  t..  recedes  from  the  current  time  t  .  In  addition, 
we  also  fade  the  uncertainly  associated  with  the  initial  state  estimate 
y(o,n). 


(y(0,n)-y(0)(y(0,n)-y(0))TJ  P(0)exp^  \  c.^ 


(73) 


With  the  above  definitions  the  fading  memory  Kalman  filter  estimates 
y(i,n)  may  be  derived.  Actually,  we  only  consider  the  estimates  y(n,n) 
and  drop  the  double  subscript.  We  will  not  derive  the  fading  memory 
estimates  here  but  will  only  indicate  the  changes  from  the  usual  Kalman 
filter  equations.  The  only  changes  occur  in  the  equations  for  computing 
the  predicted  covariance  matrices  which  become 


P(kjk-1)  =  a(k-l)A(k-l)P  (k-l)Al(k-l)  *  Q  (k-1) 

XX 


(7A) 


Pb(k|k-1)  =  a(k-l)B(k-l)Pb(k-l)BI(k-l) 


PJ{b(k|  k-1)  =  a(^-l)A(k-l)Pxb(k-l)BT(k-l) 


(75) 

(76) 


where 


a(k-i)  =  e 


'k-1 
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The  fading  a(k)  is. chosen  to  reduce  errors  in  the  bias  estimates 

caused  by  mismodeling  while  the  state  noise  covariance  Q^(k)  is  used  to 

reduce  the  effects  of  mismodeling  in  the  process  dynamics.  The  fading 

* 

also  effects  the  weighting  for  the  zero-biase  state  estimates  x  (k),  but 
the  effect  is  normally  small  since  the  fading  factors  c^  are  nearly 
unity  in  our  application  to  instrument  bias  estimation.  The  fading 
memory  filter,  although  useful  in  reducing  errors  due  to  mismodeling  in 
bias  estimation,  is  not  entirely  satisfactory  because  it  does  not  allow 
for  treating  each  bias  term  individually. 

PROBLEMS  ENCOUNTERED  IH  BIAS  ESTIMATION.  The  major  problems 
encounterd  in  estimating  measurement  error  bias  are  concerned  with 
modeling  not  only  of  the  instrument  biases  but  also  modeling  of  the 
trajectory.  Trajectory  modeling  causes  problems  in  bias  estimation 
since  if  an  adequate  representation  of  the  trajectory  dynamics  has  not 
been  included  in  the  Kalman  filter,  it  will  be  impossible  to  separate 
the  resulting  trajectory  estimation  errors  from  the  measurement  bias 
estimates.  For  example,  if  a  constant  acceleration  is  assumed  for  the 
dynamics  of  a  missile  trajectory  when  in  reality  the  missile  is  turning, 
estimates  of  accelerometer  biases  will  be  significantly  effected  by  the 
trajectory  estimation  error  caused  by  the  constant  acceleration 
assumption.  Although  the  measurement  bias  estimates  will  be  useless  in 
cases  where  they  are  severely  confused  with  the  trajectory  modeling 
errors  the  trajectory  state  estimates  obtained  by  including  the  bias  will 
often  be  a  significant  improvement  over  the  trajectory  state  estimates 
provided  by  the  zero-bias  filter. 

Another  source  of  confusion  of  the  bias  error  estimates  occurs 
from  modeling  the  random  noise  characteristics  of  the  measurements.  If 
the  measurement  noise  is  predominantly  low  frequency  and  the  filter 
assumes  the  noise  to  be  purely  random  (usually  the  case),  the  measure¬ 
ment  bias  estimate  will  include  a  significant  portion  of  the  low  frequency 
noise  components . 
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Another  type  of  modeling  problem  occurs  in  modeling  the  measurement 
biases.  As  an  example  consider  the  following  bias  error  model  for  a 
radar  azimuth  measurement. 

AA  =  aQ  +  a^anEgSinAg  +  ajtanEgCOsAg  +  a3tanEQ  +  a^secEg 

M 

*a5^0  *  a6^*0  *  a7^0 

where  A^  and  Eg  are  measured  azimuth  and  elevation  and  Aa  is  the  azimuth 
bias  error 

aQ  =  Zero  set  error 

.  a^tanEgSinAg+a^anEgCOsAg  =.  Mislevel  error 
a3tanEQ  =  Orthogonality  error 
a^secEg  =  Collimation  error 

a,.  A  =  Servo  velocity  error 

•  • 

a^A  =  Timing  error 

•• 

a^A  =  Servo  acceleration  error 

Many  more  terms  could  be  included  in  the  above  error  model  but  the 
terms  listed  will  serve  to  illustrate  the  desired  points.  The  first 
problem  with  the  azimuth  error  model  is  obvious;  the  timing  and  servo 
velocity  error  terms  have  exactly  the  same  form  which  implies  that  a,, 
and  a6  cannot  be  separately  estimated.  Thus,  we  will  only  be  able  to 
estimate  the  sum  of  a*,  and  a^.  How  suppose  that  the  trajectory  being 


estimated  is  approximately  a  level  flying  aircraft  and  that  the  ground 
range  from  the  radar  is  such  that  the  elevation  angle  is  nearly  constant 
over  large  segments  of  the  trajectory.  In  this  case  tanEg  and  secEg  are 
effectively  zero  set  error  terms.  In  this  case  there  is  little  chance 
of  obtaining  meaningful  estimates  of  ag,  a^,  and  a^.  The  three  terms  may 
be  lumped  and  a  single  zero  set  error  estimated  or  if  the  an,  a,,  and  a, 
are  separately  estimated,  the  sum  ag+a2tanEg_a^secEg  may  be  a  meaningful 
estimate.  This  example  shows  that  each  term  in  an  error  model  must  be 
well  exercised  along  the  trajectory  in  order  to  obtain  meaningful  esti¬ 
mates  of  the  error  coefficients. 

Another  problem  in  modeling  measurement  bias  ss  involves  the  large 
number  of  error  coefficient  terms  which  would  be  required  to  adequately 
model  all  measurement  bias  errors  in  some  WSMR  missions.  Consider  again 
the  LANCE  example  given  in  the  introduction  where  there  were  66  different 
measurements.  If  an  average  of  three  error  coefficients  are  required  to 
model  a  bias  error  in  this  example,  which  is  certainly  not  unrealistic, 
198  bias  state  variables  are  required  in  the  Kalman  filter.  Even  though 
the  technique  presented  in  this  report  has  greatly  speeded  up  the 
computation  of  bias  error  estimates,  it  is  doubtful  that  the  estimation 

p 

of  198  states  is  feasible. 

To  summarize  the  above  remarks  on  modeling  problems  the  following 
requirements  must  be  satisfied  in  order  to  do  a  completely  satisfactory 
job  of  bias  estimation. 

1*  An  adequate  representation  of  the  trajectory  dynamics  must  be 
Included  in  the  filter. 

2.  A  realistic  model  of  the  measurement  noise  including  serial 
correlation  must  be  available. 


3.  Each  term  in  a  measurement  bias  error  model  must  be  well  exercised 
along  the  trajectory. 

4.  The  number  of  bias  state  variables  must  be  small  enough  for  the 
computational  problem  to  be  feasible  yet  large  enough  so  that  no  severe 
aliasing  errors  are  present. 

In  any  case  careful  study  is  required  to  determine  what  terms  to 
include  in  the  bias  error  model  and  to  determine  if  the  resulting  bias 
estimatesare  actually  due  to  the  error  source  assumed  or  are  significantly 
effected  by  some  other  factors  such  as  trajectory  modeling  errors,  low 
frequency  measurement  noise,  or  missing  terms  in  the  bias  error  model. 
Although  the  modeling  problems  presented  above  pose  some  very  difficult 
problems,  they  do  not  impose  a  serious  limitation  on  the  use  of  bias 
estimation,  if  one  is  willing  to  accept  the  premise  that  the  primary 
purpose  of  bias  estimation  is  to  improve  the  trajectory  state  estimates. 

SQUARE  ROOT  IMFLEMSTATIOK  OF  BIAS  ESTIMATION .  Several  matrix  square  root 
formulations  of  the  Kalman  filter  equations  have  been  presented  in  the 
literature.  A  comprehensive  survey  of  these  methods  along  with  an  excellent 
bibliography  is  given  in  [8],  Until  recently  we  used  the  matrix  square 
root  methods  presented  in  [2]  and  [4].  Presently  we  employ  the  square  root 
filtering  methods  described  in  [3]  which  are  summarized  below.  We  have 
found  that  the  square  root  filter  provides  computational  efficiency  as  well 
as  numerical  stability  in  the  mechanization  of  the  Kalman  filter.  The 
square  root  filtering  equations  given  below  consider  only  the  processing 
of  scalar  observations.  This  is  no  restriction  since  either  the  square 
root  method  employed  is  easily  extendable  to  the  processing  of  vector 
observations  or  the  problem  of  processing  vector  observations  may  be 
reduced  to  the  processing  of  scalar  observations. 


After  a  time  update  the  predicted  covariance  matrix  P*(k)  k-1)  of 
the  zero-bias  filter  is  decomposed  as 

P*(k|  k-1)  -  L(kJk-l)D(k|k-l)LT(k|k-l) 


where,  LT(k| k-1)  is  a  lower  triangular  matrix  having  ones  along  the 
diagonal  and  D(k|  k-1)  is  a  diagonal  matrix  having  positive  diagonal 
elements.  This  triangular  decomposition  is  related  to  the  Choleski 
decomposition.  The  algorithm  for  decomposition  is  summarized  in  Appendix 
A.  For  further  details  see  (9).  For  each  scalar  observation  occuring  at 
the  new  time  t^  an  updated  triangular  decomposition  of  the  covariance  and 
an  updated  state  estimate  are  computed.  Let  x*(k)  denote  the  zero-bias 
state  estimate  after  processing  the  ith  scalar  observation  at  t^  and  let 

(i)  (i)  (i)T 
L(k)D(k)L(k) 

be  the  triangular  decomposition  of  the  covariance  matrix  after  the  ith 
measurement  update  at  tfc.  The  updated  triangular  decomposition  satisfies 


(i)  (i)  (i)T 
L(k)D(k)L(k) 


(i-1)  (i-1)  (i-l)T 


(i)  (i)  (i)T 


L(k)  D(k)  L(k)  -  C^'y^'y; 


where  the  vector  and  the  scalar  are  computed  from 


(i-1) (i-1)  (i) 
L(k)  D(k)  u* 


C*1J  -  l/(R^(k)  +  ■ 


.  (i-l)T  T 
u*  -  L(k)  Hi(k) 

(0)  (0) 

Also  L(k)  »  L(k|  k-1)  and  D(k)  ■  D(kj  k-1) . 


The  updated  state  estimate  x*(k)  is  computed  from 
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*(i> 

x*(k) 


.  (i“l) 
x*(k) 


+  C. 


(i)  (i) 

y* 


*(i-D 

^zi(k)  ’  ui(k)x  (k)  ) 


(82) 


*<°>  * 

x  (k)  -  x  (kjk-1). 


(i)  CD 

An  algorithm  is  derived  in  [3]  for  computing  the  L(k)  and  D(k)  in  (78)  given 
(i-1)  (i-1)  (i)  (i) 

L(k)  ,  D(k)  ,  C*  ,  and  y^  .  This  algorithm  is  summarized  in  Appendix 

B. 


The  triangular  decomposition 

Pfe(k)  -  Lb(k)Db(k)Lb(k) 

is  also  used  for  the  bias  filter.  We  will  consider  only  the  case  where 
the  elements  of  the  bias  vector  are  assumed  to  be  constant.  "In  addition 
to  being  easier  to  handle  numerically  the  case  of  constant  bias  is  probably 
the  most  useful  in  practice.  The  triangular  decomposition  for  the 
more  general  case  of  linear  bias  dynamics  results  In  sHghtJy  more 
numerical  work.  With  constant  bias  dynamics  the  square  root  bias  filter 


equations  at  a  time  update  are 

b(kj  k-1)  »  b(k-l) 

(83) 

^(klk-l)  -  Lfa(k-1) 

(84) 

Db(k|k-1)  »  Db(k-1). 

(85) 

For  each  scalar  measurement  update  at  the  new  time  t^  the  updated 
triangular  decomposition  satisfies 

(i)  (i)  (i)T  (i-1)  (i-1)  (i-l)T  (i).(i)  (i)T 

Lb(k)Db(k)Lb(k)  »  ^(k)  Db(k)  Lb(k)  -  Cb  yb  yb  (86) 
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(i)  (i) 


The  vector  yb 

and  the  scalar  Cb  are  computed  from 

(1) 

(i-1)  (i-1)  (i)  ’ 

yb 

"  Lb(k)  Db(k)  % 

(87) 

(i) 

(i)  (i-1)  (i) 

.  Cb 

"  l/(CjL(k)  +  ub  Db(k)  ) 

(88) 

(i-l)T 

"  Lj, (k)  si(k) 

(89) 

c^k) 

*T  (i-D  * 

*  R^(k)  +  uA  D(k)  u± 

(90) 

(0) 

(0) 

^(k)  -  L(k|  k-1) ,  Db(k)  »  D(k|  k-1) 


The  updated  bias  state  estimate  is 


*(i)  *(i-l)  <i)  (i)  *  T  (i-1) 

b(k)  -  b(k)  +  Cb  yfa  (r^k)  -  Si(k)b(k)  )  (91) 

*  *U> 

fji  (k)  "  Z±(k)  -  Hi(k)x  (k)  (92) 

(i-1) 

8i(k)  "  Gi(k)  +  Vk)T(k)  (93) 


For  each  measurement  update  at  t^  the  combining  matrix  T  is  updated  as 
(i)  (i-1)  (i)  (i) 

T(k)  «  T(k)  -  C*  y*  sj[(k)  (94) 


T(k)  -  T(kJ k-1) 

■BIAS  FILTER  REINITIALIZATION.  The  dimension  of  the  bias  state  vector  may 
change  frequently  during  the  execution  of  the  filtering  process.  This 
dimension  change  is  required  when  a  measuring  instrument  begins  to  take 
observations  after  the  initialization  of  the  filter  or  when  an  instrument 
Is  deleted  from  the  filtering  solution  because  it  has  stopped  taking 
observations,  its  bias  is  considered  too  large,  or  its  observations  arc 
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chronically  inconsistent  with  their  statistics.  The  decomposition  of 
the  optimal  filter  into  a  zero-bias  filter  and  a  bias  filter  requires 
that, the  zero-bias  state  estimate  be  orthogonal  (in  the  usual  statistical 
sense)  to  the  bias  state  estimate.  Thi3  condition  must  be  satisfied  when 

*  .  reinitializing  the  bias  filter.  In  addition  we  require  that  there  be  no 

A 

change  in  the  trajectory  state  estimate  x  and  the  remaining  bias  state 

*  estimates  due  to  reinitialization. 


The  orthogonality  and  continuity  conditions  are  automatically  met 
when  a  new  measuring  instrument  first  enters  the  B.E.T.  program  provided 
We  assume  that  its  initial  bias  estimates  are  uncorrelated  with  the 
zero-bias  state  estimate.  However,  when  a  measuring  instrument  is  deleted 
from  the  filtering  solution,  considerable  effort  is  required  in  reinitial¬ 
ization  in  order  to  meet  the  orthogonality  and  continuity  conditions. 

Let  x*,  P*  denote  the  zero-bias  state  estimate  and  its  covariance 
just  prior  to  dropping  a  measurement  from  the  filtering  solution  and  let 
x+  and  P+  be  the  same  quantities  after  dropping  the  measurement  and 
reinitializing  the  filter.  Let  b  and  b,  be  the  bias  state  estimate  before 

*"  *T* 

A 

and  after  dropping  a  measurement,  b,  is  formed  by  deleting  the  component 

A  * 

of  b  corresponding  to  the  measurement  being  dropped.  T  and  T,  are  the 

“  "  T 

combining  matrices  before  and  after.  T.  has  one  less  column  than  T_. 

P^  and  P^  are  the  bias  covariance  matrices  before  and  after  dropping  a 

+  *• 

measurement.  P^  is  formed  from  P^  by  deleting  the  row  and  column 

corresponding  to  the  measurement  being  dropped.  The  updating  equation 
for  x*  is 

x*  «  x*  +  t1(bi  -  pTb+)  (95) 

A 

where  t.  is  the  column  of  T_  being  deleted  and  b.  is  the  bias  estimate  of 
measurement  being  dropped.  The  vector  p  is  chosen  so  that  x+  and  b+  will 
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be  orthogonal.  Using  the  triangular  decompositions  P_  »  L_D_L  and 

jl  T  ft  *” 

P*  *  L.D.Lr  the  updating  equation  for  P  may  be  written  as 

T  *T  *T*  T 

L+D+L+  "  L-D_L-  +  (pb  (i,i)  ”  PTpb  p)tiri 

"■  + 


where  P.  (i,i)  is  the  diagonal  element  of  the  row  and  column  of  P, 
b  b 

which  are  being  deleted.  The  algorithm  of  Appendix  B  is  used  to  compute 

L.  and  D,  in  the  above  equation.  Similarly,  using  the  triangular  de- 
+  +  T  T 

compositions  P^  **  and  P^  =  delete  the  ith  row  and 

-  -  -  -  +  +  +  + 

column  of  both  .  and  .  Call  the  resulting  matrices  and  D'  . 

The  updating  equation  for  P^  can  then  be  written  as 


V vv  "  LbDb.Lb.  +Wi 
+  +  +  —  —  — 


where  is  the  column  deleted  from  and  d_^  is  the  diagonal  element 

deleted  from  D,  .  L,  and  D.  are  computed  using  the  algorithm  of 
b  b.  to. 

-  +  + 

Appendix  B.  Having  computed  Lb  and  the  vector  p  which  is  chosen  so  that 

+  + 

x.  and  b,  will  be  orthogonal  is  computed  by  solving  the  triangular 

*v  T 

equation 


Lb/  •  h 


for  y  and  then  solving  the  triangular  equation 


,T  n-l 

L  p  -  D  y. 

°+  + 


In  order  to  satisfy  the  conditions  that  the  state  estimate  be 
unchanged  after  reinitialization  of  the  bias  filter,  set 


f+  -  T'  +  tjp1 


(100) 


30 


where  T*  is  formed  by  deleting  the  ith  column,  t, ,  from  T  .  This  will 

A  “  A  ^ 

make  x+“  x_.  The  bias  estimates  automatically  remain  unchanged. 


r 
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APPENDIX  A 


Let  P  be  an  nxn  symmetric,  positive  definite  matrix.  The  following 
algorithm  computes  a  lower  triangular  matrix  L  with  L(i,i)  »  1  and 
a  diagonal  matrix  D  with  d.^  **  d(i,i)>0  such  that  P  *•  LDLT.  For 
further  details,  see  [9], 

dx  »  p(l,l) 

p*(k,l)  *  p(k,l)  k>l 

Jt(k,l)  -  p*lk,l)/d1 
i-1 

d  -  p(i,i)  -l  p  (i,j)f(i,j)  k>i 


-7<y-  ,4  ^ ,  «T 


;  5  *?<£»  ,,*>v  7>?r  **&***??, <fi*s 


APPENDIX  B 


Let  L  be  a  lower  unit  triangular  matrix  and  D  be  a  positive  diagonal 

T 

matrix  D  such  that  P  -  LDL  is  positive  definite.  Given  a  scalar  c 

T 

and  a  vector  x  such  that  P"  *»  P  +  cxx  is  positive  definite,  compute 
a  unit  lower  triangular  matrix  h'  and  a  positive  diagonal  matrix  D" 
such  that  P"  =  L'D'L'  .  The  following  algorithm  computes  L'  and  D" 
given  L,  D,  c,  and  x. 


(1)  CD 

,  8  C  |  X  c  X 


K 


.  .  (1)  (i)2 

°i  +  C  Xi 


x(i+1) 

J 

A'(j,i) 


,(i) 


Mj,i) 


x<*> 

+ - — 


(i+1) 

j 


i  »  1,  n 
j  *  i+1,  n 


»> 
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